Prepare functions
# function calculates all relevant models
run_models <- function(y, lvl1_x, lvl2_x, lvl2_id, data, ctrls=F){
# subset data
data = data %>%
dplyr::select(all_of(y), all_of(lvl1_x), all_of(lvl2_x), all_of(lvl2_id),
popdens, all_of(y))
data = data %>%
dplyr::rename(y = all_of(y),
lvl1_x = all_of(lvl1_x),
lvl2_x = all_of(lvl2_x),
lvl2_id = all_of(lvl2_id)
)
# configure optimization procedure
ctrl_config <- lmeControl(opt = 'optim', maxIter = 100, msMaxIter = 100)
# baseline
baseline <- lme(fixed = y ~ 1, random = ~ 1 | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# random intercept fixed slope
random_intercept <- lme(fixed = y ~ lvl1_x + lvl2_x,
random = ~ 1 | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# random intercept random slope
random_slope <- lme(fixed = y ~ lvl1_x + lvl2_x,
random = ~ lvl1_x | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# cross level interaction
interaction <- lme(fixed = y ~ lvl1_x * lvl2_x,
random = ~ lvl1_x | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# create list with results
results <- list('baseline' = baseline,
"random_intercept" = random_intercept,
"random_slope" = random_slope,
"interaction" = interaction)
if (ctrls == 'dem' | ctrls == 'prev'){
# random intercept random slope
random_slope_ctrl_dem <- lme(fixed = y ~ lvl1_x + lvl2_x + popdens,
random = ~ lvl1_x | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# cross level interaction
interaction_ctrl_main_dem <- lme(fixed = y ~ lvl1_x * lvl2_x + popdens,
random = ~ lvl1_x | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# cross level interaction
interaction_ctrl_int_dem <- lme(fixed = y ~ lvl1_x * lvl2_x + lvl1_x * popdens,
random = ~ lvl1_x | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# create list with results
results <- list('baseline' = baseline,
"random_intercept" = random_intercept,
"random_slope" = random_slope,
"interaction" = interaction,
"random_slope_ctrl_dem" = random_slope_ctrl_dem,
"interaction_ctrl_main_dem" = interaction_ctrl_main_dem,
"interaction_ctrl_int_dem" = interaction_ctrl_int_dem)
}
if (ctrls == 'prev'){
# random intercept random slope
random_slope_ctrl_prev <- lme(fixed = y ~ lvl1_x + lvl2_x + popdens + rate_day,
random = ~ lvl1_x + rate_day | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# cross level interaction
interaction_ctrl_main_prev <- lme(fixed = y ~ lvl1_x * lvl2_x + popdens + rate_day,
random = ~ lvl1_x | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# cross level interaction
interaction_ctrl_int_prev<- lme(fixed = y ~ lvl1_x * lvl2_x + lvl1_x * popdens + rate_day,
random = ~ lvl1_x + rate_day | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# create list with results
results <- list('baseline' = baseline,
"random_intercept" = random_intercept,
"random_slope" = random_slope,
"interaction" = interaction,
"random_slope_ctrl_dem" = random_slope_ctrl_dem,
"interaction_ctrl_main_dem" = interaction_ctrl_main_dem,
"interaction_ctrl_int_dem" = interaction_ctrl_int_dem,
"random_slope_ctrl_prev" = random_slope_ctrl_prev,
"interaction_ctrl_main_prev" = interaction_ctrl_main_prev,
"interaction_ctrl_int_prev" = interaction_ctrl_int_prev)
}
if(ctrls == 'exp'){
# random intercept random slope
random_slope_exp <- lme(fixed = y ~ (lvl1_x + I(lvl1_x^2)) + lvl2_x,
random = ~ (lvl1_x + I(lvl1_x^2)) | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# cross level interaction
interaction_exp <- lme(fixed = y ~ (lvl1_x + I(lvl1_x^2)) * lvl2_x,
random = ~ (lvl1_x + I(lvl1_x^2)) | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# create list with results
results <- list('baseline' = baseline,
"random_intercept" = random_intercept,
"random_slope" = random_slope,
"interaction" = interaction,
"random_slope_exp" = random_slope_exp,
"interaction_exp" = interaction_exp)
}
return(results)
}
# extracts table with coefficients and tests statistics
extract_results <- function(models) {
models_summary <- models %>%
map(summary) %>%
map("tTable") %>%
map(as.data.frame) %>%
map(round, 10)
# %>% map(~ .[str_detect(rownames(.), 'Inter|lvl'),])
return(models_summary)
}
# calculates comparison of all models in model list
compare_models <- function(models) {
mdl_names <- models %>% names()
str = ''
for (i in mdl_names){
mdl_str <- paste('models$', i, sep = '')
if(i == 'baseline'){
str <- mdl_str
}else{
str <- paste(str, mdl_str, sep=', ')
}
}
anova_str <- paste0('anova(', str, ')')
mdl_comp <- eval(parse(text=anova_str))
rownames(mdl_comp) = mdl_names
return(mdl_comp)
}
Rescale Data
df_uk_scaled <- df_uk %>% dplyr::select(ut_area, time, pers_o,
pers_c, pers_e, pers_e, pers_a, pers_n,
popdens, rate_day) %>%
mutate_at(vars(-ut_area, -time), scale)
df_uk_scaled %>% head()
df_uk_socdist_scaled <- df_uk_socdist %>% dplyr::select(nuts3, time, pers_o,
pers_c, pers_e, pers_e, pers_a, pers_n,
popdens, socdist_tiles) %>%
mutate_at(vars(-nuts3, -time), scale)
df_uk_socdist_scaled %>% head()
NA
Predict prevalence
prevalence ~ openness
models_o_covid <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_o',
lvl2_id = 'ut_area',
data = df_uk_scaled,
ctrls = 'dem')
extract_results(models_o_covid)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_ctrl_dem
$interaction_ctrl_main_dem
$interaction_ctrl_int_dem
compare_models(models_o_covid)
NA
prevalence ~ conscientiousness
models_c_covid <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_c',
lvl2_id = 'ut_area',
data = df_uk_scaled,
ctrls = 'dem')
extract_results(models_c_covid)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_ctrl_dem
$interaction_ctrl_main_dem
$interaction_ctrl_int_dem
compare_models(models_c_covid)
NA
NA
prevalence ~ agreeableness
models_a_covid <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_a',
lvl2_id = 'ut_area',
data = df_uk_scaled,
ctrls = 'dem')
extract_results(models_a_covid)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_ctrl_dem
$interaction_ctrl_main_dem
$interaction_ctrl_int_dem
compare_models(models_a_covid)
NA
NA
prevalence ~ neuroticism
models_n_covid <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_n',
lvl2_id = 'ut_area',
data = df_uk_scaled,
ctrls = 'dem')
extract_results(models_n_covid)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_ctrl_dem
$interaction_ctrl_main_dem
$interaction_ctrl_int_dem
compare_models(models_n_covid)
NA
NA
Predict social distancing
social distancing ~ openness
models_o_socdist <-run_models(y = 'socdist_tiles',
lvl1_x = 'time',
lvl2_x = 'pers_o',
lvl2_id = 'nuts3',
data = df_uk_socdist_scaled,
ctrls = 'dem')
extract_results(models_o_socdist)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_ctrl_dem
$interaction_ctrl_main_dem
$interaction_ctrl_int_dem
compare_models(models_o_socdist)
NA
social distancing ~ conscientiousness
models_c_socdist <-run_models(y = 'socdist_tiles',
lvl1_x = 'time',
lvl2_x = 'pers_c',
lvl2_id = 'nuts3',
data = df_uk_socdist_scaled,
ctrls = 'dem')
extract_results(models_c_socdist)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_ctrl_dem
$interaction_ctrl_main_dem
$interaction_ctrl_int_dem
compare_models(models_c_socdist)
NA
NA
social distancing ~ agreeableness
models_a_socdist <-run_models(y = 'socdist_tiles',
lvl1_x = 'time',
lvl2_x = 'pers_a',
lvl2_id = 'nuts3',
data = df_uk_socdist_scaled,
ctrls = 'dem')
extract_results(models_a_socdist)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_ctrl_dem
$interaction_ctrl_main_dem
$interaction_ctrl_int_dem
compare_models(models_a_socdist)
NA
NA
social distancing ~ neuroticism
models_n_socdist <-run_models(y = 'socdist_tiles',
lvl1_x = 'time',
lvl2_x = 'pers_n',
lvl2_id = 'nuts3',
data = df_uk_socdist_scaled,
ctrls = 'dem')
extract_results(models_n_socdist)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_ctrl_dem
$interaction_ctrl_main_dem
$interaction_ctrl_int_dem
compare_models(models_n_socdist)
NA
NA
prevalence ~ conscientiousness
models_c_covid_exp <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_c',
lvl2_id = 'ut_area',
data = df_uk_scaled,
ctrls = 'exp')
extract_results(models_c_covid_exp)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_exp
$interaction_exp
compare_models(models_c_covid_exp)
NA
prevalence ~ agreeableness
models_a_covid_exp <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_a',
lvl2_id = 'ut_area',
data = df_uk_scaled,
ctrls = 'exp')
extract_results(models_a_covid_exp)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_exp
$interaction_exp
compare_models(models_a_covid_exp)
NA
prevalence ~ neuroticism
models_n_covid_exp <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_n',
lvl2_id = 'ut_area',
data = df_uk_scaled,
ctrls = 'exp')
extract_results(models_n_covid_exp)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_exp
$interaction_exp
compare_models(models_n_covid_exp)
NA
Create overview table
Define function to create overview tables
summary_table <- function(models, dv_name){
temp_df_ctrl_main <- NULL
temp_df_ctrl_int <- NULL
for (i in models){
results <- i %>% extract_results()
results_ctrl_main <- results$interaction_ctrl_main_dem['lvl1_x:lvl2_x',]
temp_df_ctrl_main <- temp_df_ctrl_main %>% rbind(results_ctrl_main)
results_ctrl_int <- results$interaction_ctrl_int_dem['lvl1_x:lvl2_x',]
temp_df_ctrl_int <- temp_df_ctrl_int %>% rbind(results_ctrl_int)
}
names_ctrl_main <- paste0(dv_name, '~', c('o', 'c', 'e', 'a', 'n'), '*time', '_crtl_popdens')
names_ctrl_int <- paste0(dv_name, '~', c('o', 'c', 'e', 'a', 'n'), '*time', '_crtl_popdens*time')
rownames(temp_df_ctrl_main) <- names_ctrl_main
rownames(temp_df_ctrl_int) <- names_ctrl_int
sum_tab <- rbind(temp_df_ctrl_main, temp_df_ctrl_int) %>% round(4)
return(sum_tab)
}
Create overview tables
# prevalence
models_prev <- list(models_o_covid,
models_c_covid,
models_e_covid,
models_a_covid,
models_n_covid)
sum_tab_prev <- summary_table(models_prev, dv_name = 'prev')
write.table(sum_tab_prev, quote=F)
Value Std.Error DF t-value p-value
prev~o*time_crtl_popdens 0.0037 0.0038 2413 0.9907 0.3219
prev~c*time_crtl_popdens -0.0177 0.0037 2413 -4.7421 0
prev~e*time_crtl_popdens -0.0033 0.0038 2413 -0.8744 0.382
prev~a*time_crtl_popdens -0.0131 0.0037 2413 -3.5102 5e-04
prev~n*time_crtl_popdens 0.0138 0.0037 2413 3.6925 2e-04
prev~o*time_crtl_popdens*time -0.0026 0.0039 2412 -0.6683 0.504
prev~c*time_crtl_popdens*time -0.0042 0.0051 2412 -0.8169 0.4141
prev~e*time_crtl_popdens*time -0.0049 0.0037 2412 -1.321 0.1866
prev~a*time_crtl_popdens*time -0.0048 0.0041 2412 -1.1739 0.2406
prev~n*time_crtl_popdens*time 0.0074 0.0039 2412 1.8814 0.06
# social distancing
models_socdist <- list(models_o_socdist,
models_c_socdist,
models_e_socdist,
models_a_socdist,
models_n_socdist)
sum_tab_socdist <- summary_table(models_socdist, dv_name = 'socdist')
write.table(sum_tab_socdist, quote=F)
Value Std.Error DF t-value p-value
socdist~o*time_crtl_popdens -4e-04 0.0018 2396 -0.239 0.8111
socdist~c*time_crtl_popdens 8e-04 0.0018 2396 0.4223 0.6728
socdist~e*time_crtl_popdens -0.0012 0.0018 2396 -0.7034 0.4819
socdist~a*time_crtl_popdens 2e-04 0.0018 2396 0.1142 0.9091
socdist~n*time_crtl_popdens 2e-04 0.0018 2396 0.1124 0.9105
socdist~o*time_crtl_popdens*time 0 0.0019 2395 0.0154 0.9877
socdist~c*time_crtl_popdens*time -7e-04 0.0024 2395 -0.2708 0.7866
socdist~e*time_crtl_popdens*time -0.0011 0.0018 2395 -0.6193 0.5358
socdist~a*time_crtl_popdens*time -6e-04 0.002 2395 -0.3054 0.7601
socdist~n*time_crtl_popdens*time 9e-04 0.0019 2395 0.4509 0.6521
Social distancing